MICRO-MACRO MODELLING OF AN ARRAY OF SPHERES 
INTERACTING THROUGH LUBRICATION FORCES 
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Abstract. We consider here a discrete system of spheres interacting through a lubrica- 
tion force. This force is dissipative, and singular near contact: it behaves like the recip- 
rocal of interparticle distance. We propose a macroscopic constitutive equation which is 
built as the natural continuous counterpart of this microscopic lubrication model. This 
model, which is of the newtonian type, relies on an elongational viscosity, which is pro- 
portional to the reciprocal of the local fluid fraction. We then establish the convergence 
in a weak sense of solutions to the discrete problem towards the solution to the partial 
differential equation which we identified as the macroscopic constitutive equation. 

Resume. Nous considerons ici un systeme discret de spheres en interaction a travers 
une force de lubrification. Cette force est dissipative et singuliere pres du contact : 
elle se comporte comme l'inverse de la distance interparticulaire. Nous proposons une 
equation constitutive macroscopique qui est construite comme le pendant continu de ce 
modele discret de lubrification. Ce modele, de type Newtonien, repose sur une viscosite 
elongationnelle proportionelle a l'inverse de la fraction locale de fluide. Nous etablissons 
ensuite la convergence dans un sens faible des solutions du probleme discret vers les 
solutions de l'equation aux derivees partielles que nous avons identifiee comme l'equation 
macroscopique constitutive. 
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1. Introduction 



Description of the macroscopic behaviour of diluted suspensions of rigid particles in a 
Newtonian fluid has given rise to a large amount of papers over the last century. The first 
models (see [6], [81 §22] where Eintein's approach is detailed, or [5]) were based on the 
asumption that particles do not interact, which restricts their domain of validity to highly 
diluted suspension. This approach was extended to semi-diluted suspension ([!]), which 
leads to second order asymptotic expansions of the apparent viscosity with respect to the 
solid fraction. 

More recently, some authors have investigated the other end of the dilution scale: the case 
of highly packed suspensions of rigid spheres. Direct or semi-direct numerical simulations 
have been carried out to investigate the behaviour of highly concentrated fluid-particle 
mixtures ([2], [9]). 

In |10j . a first attempt was proposed to investigate the behaviour of the apparent shear 
viscosity of a suspension in the neighbourhood of the maximal packing solid fraction $rnax* 
A model is proposed, which gives a shear viscosity which is of the order (1 — f/$ max ) _2 
where $ is the solid fraction. In [4], the authors investigate the asymptotic behaviour, as e 
goes to 0, of a set of particles under the assumption that distances between neighbouring 
particles is subject to behave like e. In this framework, the authors establish that the 
apparent shear viscosity behaves like l/e 3//2 . This approach extents a previous work ([7]) 
where periodic arrays of spheres are considered. In this context, the elongational viscosity 
can be shown to behave like 1/d where d is the constant distance between neighbouring 
spheres. 

The approach we propose here is based on a simpler model from the geometric standpoint, 
as the spheres are supposed to be aligned. On the other hand it generalizes these works in 
the sense that no assumption is made on the distances: the macroscopic behaviour depends 
on the solid fraction only. Contacts between neighbouring particles are even allowed, and 
a special attention has been paid to the way we express the continuous model so that 
macroscopic clusters can be taken into account (the local viscosity within a cluster is 
infinite). 

We prove, as expected, that the limit elongational viscosity behaves singularly with re- 
spects to the vanishing fluid fraction. This approach leads to an equation of the elliptic 
type 



where v is the velocity field, p the solid fraction (which is 1 when all particles are in 
contact), and / an external body force. 




2. Discrete model 



Consider two rigid spheres imbedded in a viscous fluid, subject to move horizontally. 
Denoting by q\ and q2 the abscisses of their centers, by u\ and U2 their instantaneous 
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Figure 1. Lubrication model 



velocities and by d the border-to-border distance, the leading term in the asymptotic 
expansion of the interaction force is (see [3]) 



f1\ TP ui-ui 

(1) *U 2 = -k— 5-, 

where k is a constant which depends on the viscosity of the lubricating fluid and the radii 
of the spheres. We shall take k = 1 in what follows. Consider now an array of N + 1 
spheres, on the x-axis, with the same radius e. We set the first and the last sphere at 
position and 1, respectively. As a consequence, the number of degrees of freedom is 
N — 1, whereas the number of actual spheres is N + 1. 



go qi qt-i q% qN-i qN 




Figure 2. Geometry 

Definition 2.1. Given a vector of positions q = (qi)i<i<N~i, we say that q is e-feasible 
(spheres do not overlap), if 

qi > 2e , QN—1 < 1 - 2e , % - - 2e > V* = 1, N, 

and strictly e-feasible if all inequalities are strict (spheres do not touch). 

We denote by di = qi — — 2e the distance between spheres i and i — 1, by u, the 
instantaneous velocity of sphere i, and by u = (uj)i<j<jv_i the velocity vector. Velocities 
of the extremal spheres and N + 1 are taken as (see remark 12.41 for non-zero extremal 
velocities). Given a strictly e-feasible vector q, we define A(q) as the (N — 1) x (N — 1) 
tridiagonal stiffness matrix 
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with di = qi — qi-\ — 2e. Consider now a set of forces /i, /2,- • • /iV-ii an d the corresponding 
vector f. From ([!]), the balance of forces reads 



(3) -A(q)u + f = 0. 

Proposition 2.2. Given a strictly e-feasible vector q € IR^ -1 , a force field f € IR^ -1 , 
problem {3p /ias a unique solution u, and we shall write 

u = («i)i<i<jv-i = ^(q,f,£)- 

T/iis solution can be expressed 

if i 1 

(4) = — UD N -Di)J2D k f k + Di £ (D N -D k )f k \ Vi = l...JV-l 

j 



fc=i fc=t+i 



Proof. Matrix A, which is similar to the matrix obtained by discretizing the Laplace 
operator with Dirichlet boundary condition by finite differences, is symmetric positive 
definite, and the vector u is immediately checked to solve the system. □ 



This approach can be extended to e-feasible situations in a large sense (particles are 
allowed to get into contact). As the interaction force (which tends to penalize the relative 
velocity) blows up when particles tend to get into contact, we simply consider that two 
particles in contact have the same velocity. This situation can be formalized the following 
way: The N + 1 particules form N c clusters, and the k-th cluster contains the N k + 1 
particles i k , i k + 1, . . . , i k + N k (see Fig. [3]). The balance of forces now reads 



(5) Vi<£u k [i k ,i k + N k ], 

(6) Vke[l,N e ], { 



U i+ i -Ui Ui — Ui-1 



di+l di 

■ ■ = U ik+N k , 
Ui k +N k +1 ~ U ik+Nk U ik - U ik ^i 



Ui k — Ui k+1 



d 



i k +N k +l 



di k 



i k +N k 
i=i k 



2 s 



Qi k -1 Qi k H k +N k <iik+N k +l 

Figure 3. Non-strictly e-feasible configuration 

Proposition 2.3. Given an e-feasible vector q £ IR^ -1 , a force field f G IR^ -1 , prob- 
lem {5p(0|) has a unique solution, and we shall write as before 

u = (ui)i<i<Ar_i = 7>(q,f,e). 
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An explicit expression of this solution is given by Q). 

Remark 2.4. It is possible to set extremal velocities uq and it/v to non-zero values : in 
that case, the balance of forces is given by 

A(q)u = f + b, 

where b contains the non homogeneous dirichlet conditions : b = (ug/di, 0, . . . , 0, ujv/d/v)*- 
The extension of Proposition 12.31 to that case is straightforward. 



3. Micro to macro modelling: a heuristic approach 



The purpose of this section is to derive heuristically a macroscopic constitutive equation 
from the discrete balance of forces. 

We consider an e- feasible configuration of N+l particles (see definition ^. 1|) q. We suppose 
that uq and un are given, and the external force is taken 0. The balance of forces reads 

(7) A(q)u = b, 

where A is given by ([2]) and b contains the non homogeneous dirichlet conditions b = 
(u /dx,Q, . . . ,0, wjv/djv)*- 



The force exerted on the Nth particle by the others is 

J A 



f sys _ UN-1 — UN 



As the solution to ([7]) can be expressed 



Ui = u + (itjv - "o)i 
L) N 



it comes 



(8) f N =-^~, 

JV 

where Djy = dj is the quantity of vacuum between qo and qw- 

3=1 

Given this discrete stress tensor, we shall conjecture the continuous one. Suppose the 
density of particles is high and note the vacuum fraction in the neighbourhood of a point 
x by D(x). The force exerted on point x by the part of the system on its right side is 
given by 

F { right side of x}-+{x} = ^ F [ x ' x +v}^{^}- 
Using ([8j), it can be expressed by 

_ V(X + T])-V(X) 

r {right side of x}-+{x} ~ ^"o jj^dx ' 
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which gives the following continuous stress tensor: 

F {right side of *}-{*} ~ D(x) dxV ^' 

Following a very standard approach, we establish now the resulting equation which ex- 
presses the local balances of forces overall a one-dimensional domain. We denote the solid 
fraction by p and suppose that an external force with density / is exerted on the system. 
We express the balance of forces on a part f2 = [a, b] of the system. The external force 
exerted on f2 is given by 

Fext = I Pf 

Jn 

and the force exerted by the rest of the system on SI is 



_ d x v{b) d x v(a) _ f / 1 
Kys ~ D{b) D(a) ~ L x \D( x f xV[x) 



so that the balance of forces on f2 reads 



This last result being true for all O we have 

-d x (jj d x v ! _ /'/ 
and, since D = 1 — p, we finally obtain 

(9) - d x {rzTp 9 ^ ) Pf- 



4. Asymptotic behaviour of the discrete solutions 



Let / denote ]0, 1[. Firstly, we build a new operator V, which is our key tool to connect 
the microscopic level to the macroscopic one. This operator is defined the following way: 
Given e > 0, an e-feasible position vector q (as stated by Definition 12. 1\ it represents a 
distribution of particles whose centers are located at i = 1, . . . , JV, with common radius 
e), a force density / € L 1 ^), we define vector u as "P(q, f £ ,e) (see Proposition 12.21 or 12.31 
depending on whether q is strictly feasible or feasible), where f £ is defined by 

ft = Ij- f ^ ds for 1 < i < iV - 1. 

Now, to this vector u = "P(q, f £ ,e) £ IR^ -1 , we associate a piecewise affine function u 
defined by 

(10) u G C°(J) , u affine on [q t , q i+x ] Vi = 0, . . . , JV - 1 , u( qi ) = m V* = 0, . . . , JV. 



We shall write it = "P(q, /, e). 
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In the same spirit, for any e > 0, and any e-feasible position vector q G R N 1 , we define 
x(q, e) as the characteristic function of the solid phase associated to the (q, e) distribution: 



(11) 



N-l 

X(q,e) = V 1 



7 j 

8=1 



+ l]0,e[ + ijl-e,!!- 



2 s 



Figure 4. Definition of ?/ 



Before we state the main convergence theorem, we still have to give a sense to Q when 
the density p G [0, 1] is allowed to take value 1, even on a set of non-zero measure. 

Proposition 4.1. Let K : 1 1— > IR + U{+oo} 6e measurable, K(x) > a > /or almost every 
x G J ; 99 € H~ l (I), and let J be defined as 

v G FqV) ' — ► J(w) = J K(x)\d x v\ 2 - < <p,v > G [R + U{+oo}. 

There exists a unique u G Hq(I) which realizes the minimum of J over Hq(I). If there 
exists f G L 1 suc/i £/iaf < tp,v >= J fv, we shall say that u is a generalized solution to 

-d x {K(x)d x u) = f. 



Proof. The functional J is convex (strictly convex over its domain), coercive, and it can 
be written 

J{v) = sup ( / mm(K (x) , n)\d x v\ 2 — / fv) , 
n€N \Jl J I ) 

thus it is l.s.c. as a supremum of a family of l.s.c functions. Therefore it admits a unique 
minimizer. Note that the minimization problem is equivalent to the problem which consists 
in minimizing the same functional J over the set 



(12) H K = \ v G Hl(I),d x v = a.e. in D(K) C and / 

{ Jd{k) 



K(x)\d x v\ 2 < 00 } . 
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where D(K) = {x € I, K{x) < +00} is the domain of K. Consequently, u is characterized 
by the variational formulation 



(13) 



u € Hk , I K(x)d x ud x v 

'D(K) 



fv Vu € Hk- 



□ 



We may now state the convergence result. 

Theorem 4.2. Let f € 1^(1) be a force density, and p E L°°{I) a solid fraction, with 
p{x) € [0,1] a.e. in I. Let (q e ) e be a sequence of e- feasible position vectors (see Defini- 
tion \2.1\) , q e G ¥r e with N £ = l/e. We introduce x £ = xitf i £ ) ( see $H\))> an d we assume 
that x £ converges toward p in L°°(L) weak-*. 

Then u £ = V((i £ ,f, e) solution to the discrete model (see U0\) ) converges weakly toward u 
in Hq(L) as e converges to 0, where u is the solution to 



(14) 



d x 



1-p 



d x u = pf, 



in the sense of Proposition \4^1\ (i.e. characterized by A13\) ). 



Proof. The proof is based on some technical lemmas. For readability reasons, we postpone 
the proofs of the lemmas to the end of the section. 



As a first step, we define p e , piecewise constant, as the proportion of solid on each subin- 
terval [<3f_i,<jf] in the following way : let di be the distance between particles i — 1 and i, 
then 



Vi = 0, . . . , iV e - 1 , p £ = p\ = 1 



a i+l 



on [gf_i,gf]. 



Pi k +N k +1 



1i k +N k +l 



Figure 5. p e 
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Define now w £ , constant by part, as the discrete counterpart of d x u/(l — p) given on 
[«f_i,«i] by 

us) I f-rh 



■d x u e , if pf < 1 



w e = wf = /3| , if pf = 1 

where (Pi)i k <i<i k +N k corresponding to the kth cluster, is the solution to the following 
system: 



(16) 



i k +l 



Pf 



i k +2 



P: 



ik+N k 
U ik+N k +1 ~ U i k +N k 
a i k +N k +l 



Pf 



i fc +l 



Pf 



ik+N k -l 

>£ 

ik+N k 



PI 



-^ftk^ 



-2e/ e 



ifc+JVfc-l' 
ik+Nk' 




1i k +N k 



1i k +N„+l 



Figure 6. w e 



Note that, summing up all equations of (|16p we recognize the balance of forces on the kth 
cluster given by (J6]). 

Remark 4.3. The idea behind the above construction is that j3{ can be seen as the 
cohesion force between particles i — 1 and i. A first way to notice it is to note that u £ is 
the limit of u £ ' v where u £,r? is the solution to system ([3]) with d% = rj > for i between 
ik + 1 and + iVfc and that we have 

e,rj 



Vfc, Vj€[l,iV fc ], /3f fc+j = lim 



£,77 



Another way to define these cohesion forces is to consider the following minimization 
problem : minimize 



Jty) 



i£U k [i k +l,i k +N k ] 



1 (Vj-Vj-i) f £ 

2 df ^ J * 



V, 



i=l 
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over K = {v , Vi £ ^ki^k + Ij^fc + Nk], Vi = Vi-i}. This problem is equivalent to 
problem ©([6]) and (3 £ turns up as the Lagrange multiplier associated to the constraint 

The next lemma shows that (|14p is true at the discrete level : 
Lemma 4.4. For any e > 0, 

(17) - d x (w £ ) = 2ef 

N e -1 

in H^ 1 (F), where f £ = ff$qf f<% ^ s the Dirac measure at point q £ ). 

i=i 

The idea of the proof is now to let e go to in (fTT|) . To that purpose, we first study p £ 
and 2ef £ when e tends to zero: 

Lemma 4.5. p £ — p in L°°(I). 

Lemma 4.6. 2ef £ — pf in H^ 1 ^). 

Let us now establish that u £ is bounded in Hq(I). Applying Lemma l4.4l to the test function 
u £ gives, using that q £ is e-feasible, 

11^11^1(7) < ||« £ |k°°(/)||/||Li(/) < Cll^ll^i^H/llii^). 

Hence, {u £ ) £ is bounded in Hq(I) and we can extract a subsequence of (u £ ) £ (still denoted 
by (u £ ) £ ) such that u £ — 1 u in Hq(I). In order to pass to the limit in the left-hand side 
of (fT7|) . we are going to prove that there exists a subsequence of (w £ ) e converging to a to 
in L l (I). This will follow from the next lemmas 

Lemma 4.7. (w £ ) £ is bounded in L°°\I). 

Lemma 4.8. (w £ ) £ has uniform bounded variation. 

From these two lemmas it follows that (w £ ) £ is bounded in the space of functions of 
bounded variation BV(I) and by compact injection of BV(I) in L 1 (/) we can find w £ 
L l (I) n BV{I) and a subsequence of (w £ ) £ (still denoted by (w £ ) £ ) such that 

(18) w £ — > w in L}(I) and a.e. 

Then, convergence of w £ and p £ make it possible to establish the following lemma: 
Lemma 4.9. 

d x u £ ->> (1 - p)w in L 2 {I). 

We now come to the last step of the proof : let e tend to zero in (|17p and obtain asymp- 
totically 

-d x ^JZ~ p d ^ = Pf> 
in the sense of Proposition 14.11 characterized by (|13p 
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First, Lemma I3~9l implies d x u = (1 — p)w, so u G #i/(i-p) (defined by (fl2|) ). Then, by 
Lemma 14.4 



y = (2e/ e ,v) W G iJ x (/). 
Passing to the limit on e gives, using (|18p . Lemma 14.71 and Lemma 14.6 

J wv' = J pfv V«6lf^(I). 



Finally, for any u G if 1 /( 1 _ p ), 



7/ J/ Vp^l 



(1 - p)w , f d x u , 
-v = - V 



Ip^l 1 - p Jp^l 1 - p 

and we conclude that u is the solution to (EG 



So, we proved that there exists a subsequence of (u e ) e converging to u as e tends to zero. 
Since the same work can be done for each subsequence of (u £ ) £ , we conclude that (u £ ) £ 
itself converges to u, which completes the proof of the theorem. □ 

Proof of Lemma [4.41 —d x (w £ ) = 2ef £ in H^ 1 ^). 

First, a simple computation shows that 

/,,£ _ 7 ,£ v £ -v £ \ N c i k +N k -l 

w= E {^^--^)^ + Y. E m-eu)^ 

i £ l..N £ - l 1+1 1 k=1 i=i k+ 1 

di > 0, d i+1 > 

, V/flE ^ffc ~ U f fc -! \ . ^ f uf +N +1 ~ Uf +N \ 

+ h v ftt+1 — ~) ,s * £ v <w.« — ft " +A "J 

in the sense of distributions. Then, combining this with system ([5]) ([6]) and with the 
definition of (fl6l) . we get 

^ £ = - E 2e ^<%- 
i=i 

By density of Cq°(7) in Hq(I) and continuous injection of Hq(I) in C (I), this result holds 
in H~ l {I) as required. 

Proof of Lemma 14.51 : p £ — ^ /) in L°° (I) . 

Since p e — p = (p £ — x £ ) + (x e ~~ P) an d X £ ~ ~^ P m the result will follow provided 

we show that 



x e ¥> - ffip) =0. 



1 
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By density of the set of stairs functions in and by using the fact that (p e ) e and 

{x £ )e are bounded in L°°(I), this in turn will follow from 



V(p piecewise constant on / , lim ylX l P~IP i P) = Q- 
To show this, it suffices to prove that 

r/3 



Ma, (5 < a < j3 < 1 lim IJ X s ~ j P £ j = 0. 

In order to do so take a and (3 such that < a < (5 < 1 and denote the particules whose 
center is in [a, /?] by qf , qf Q+1 , . . . , q £ jQ . Since j q t\x £ ~ P £ ) = 2e - (qf +1 - qf - df +1 ) = 
for 1 < i < N e — 1, we have 

x e - \ p £ = (x £ -p £ )~ / (x £ -p £ ). 



Then, a simple computation shows that the left-hand side converges to zero as e tends to 
zero wich completes the proof of the lemma. 

Proof of Lemma [Ml : 2ef e pf in 

By injection of Hq(I) in C°(I) it suffices to show that 

Vy> G C°(J) , lim ^2e £ /? <p\ -^pf<P = 0- 
Moreover, we can write 

JV £ -1 /JV £ -1 \ 

2£ E " ^ = E -X £ f)+ (X £ f ~ Pf) ■ 

i=l \ i=l / 

Then, using the fact that x e Z 5 ) the required result will follow as soon as we prove that 

V99 E C°(J) , limA £ = 0, 

where 

AT S -1 



To obtain this, merely compute 

We— 1 rqf+e re rl 

A £ = E / -f]f- <Pf~ <Pf 

i=1 Jqf-e Jo Jl-e 

and use uniform continuity of 92 together with the fact that q e is e-feasible to define e(5) 
by 

e{5) = sup \<p(x) -<p{y)\ 

x,yel,\x-y\<8 
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and show 

\A £ \<e(e)f2 \f\+ f £ Mf\+ t MI/I <e(e) 11/11^/)+ f W\\f\+ f Ml/1- 

i=i J ii- e J° Jl - £ J° Jl ~ e 

The result follows from this by letting e — > and using the fact that ipf € L}(I). □ 
Proof of Lemma 14.71 {w e ) e is bounded in L°°(I). 

By (fT5j) . ui e is piecewise contant, equal to wf on [gf_ l5 qf], where = y^t5 x u £ = — J^'" 1 
if pf < 1, and wf = (3f otherwise. Note that, by fp76|) and the fact that q e is e-feasible 



Vk = l...N c , Vj = l...N k , \(3f k+j \ 



ik+j-l 



< K-iI + II/IIlu/)- 



Therefore, to show that ||i,oo(j) is bounded it suffices to obtain an upper-bound for 
( w t)i s.t. pi<i ■ I n or der to do so, a simple computation gives using (j3J) 



2V e -l 



k=i+l ~ k=l 

Then, from — % — — < 1, -J^- < 1 and the fact that q e is e-feasible it follows that 



N e -1 

H\ < ]T < 11/11^(7) 

k=l 

and we conclude that 

ikiil-co < mhni), 

which completes the proof of the lemma. □ 
Proof of Lemma 14. 8t (w e ) e has uniform bounded variation. 
We recall that the variation of a function is defined by 

Var(/) = sup Ujip' : <p G C X C {I), \\<p\\ L oo {I) < l| . 

By Lemma [Q] applied to any test function tp € C^(I) and the fact that q e is e-feasible, 
we see that 

(19) Var(u, £ ) < 

so w £ has uniform bounded variation. □ 

Proof of Lemma 14.91 : d x u e — (1 — p)w in L 2 (I). 

Writing 

8 x u e - (1 - p)w = (1 - p £ )w e - (1 - p)w = {((1 - p e ) - (1 - p)) w} + {(1 - p')(vf - w)} , 
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we shall prove the weak convergence to zero in L 2 (I) of both term of the left-hand side. 
For the first term, take ip in L 2 {I). Using that w € BV(I) C L°°(I), it follows that 
wip £ 1^(1) and by Lemma 14.51 

lim / ((l-ff)-(l-p))w<p = Q, 

as required. 

We shall now prove the convergence of the second term to zero. By density of Cg°(I) 
in L 2 (I) and Lemma 14.71 it suffices to take (p in Cq°(J). Therefore, the result follows 
immediately from 



together with ((TH 



w)ip 



< 



\w 



w\ 
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